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METHOD FOR ANALYZING MM DIFFUSION DATA 
Background of the Invention 

5 

Reference to Related Applications 

This application claims priority to U.S. Provisional Application No. 
60/282,033, whose disclosure is hereby incorporated by reference in its entirety 
into the present disclosure. 

1 0 Field of the Invention 

The instant invention is directed generally to the field of magnetic 
resonance imaging {MRl), image display, image processing and image filtering 
to enhance structure and reduce noise. The present invention relates 
specifically to diffusion-weighted magnetic resonance imaging, more 
1 5 particularly to a quantitative and statistically robust numerical method for 
measuring anisotropic diffusion in an object, and most particularly to 
characterizing multiple fiber voxel data. 

Description of Related Art 

The background will be described with reference to related art 
20 publications listed in the attached appendix. All references cited are 
incorporated herein by reference. 

Magnetic resonance imaging 

Magnetic resonance imaging (MRl) is one of several approaches used 
to image physiological structures. For example, MRl has revolutionized 
25 radiological imaging of the internal structures of the human body, because of 
its noninvasiveftess and no known health hazards attributed to its use. 
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Basically, MRI involves the exposure of tissue to a variety of different 
polarizing magnetic and radio-frequency electromagnetic fields. The tissue's 
atomic nuclei respond to the fields and are subsequently processed to produce 
an image of the specimen. Nuclei exhibiting magnetic moments (spins) align 
5 themselves with the polarizing field. Depending upon the field's strength and 
the magnetogyric constant of the specific nuclear species involved, the nuclei 
precess about the polarizing field at an angular frequency (Larmor frequency). 

The spins exhibit a net magnetic moment in the direction of the 
polarizing field even though the magnetic components of the spins cancel each 

1 0 other in a plane perpendicular to the polarizing field. The net magnetic 
moment can be tilted by applying an excitation field perpendicular to the 
polarizing field and at a frequency near the Larmor frequency. The tilted 
magnetic moment comprises a transverse component rotating at the Larmor 
frequency in the plane perpendicular to the polarizing field. The magnitude 

15 and duration of the excitation field determines the extent to which the magnetic 
moment is tilted and the magnitude of the net transverse magnetic moment. 

Once the excitation field is removed, an external return coil senses the 
field associated with the transverse magnetic moment, producing a sinusoidal 
output at the Larmor frequency and an amplitude proportional to that of the 
20 transverse magnetic moment. The net magnetic moment gradually reorients 
itself with the polarizing field when the excitation field is removed resulting in 
the amplitude of the return coil output decaying exponentially with time. 

The rate at which the return coil output decays is dependent upon, and 
indicative of, the composition of the specimen. If the excitation field has a 
25 broad frequency band, the return coil output may include components 

associated with the transverse magnetic components of a greater variety of 
frequencies. A Fourier analysis of the output allows the different frequencies, 
indicative of different chemical or biological environments, to be distinguished. 
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Employing an excitation field that has a narrow frequency band, excites only a 
relatively narrow band within a nuclear species, resulting in the transverse 
magnetic component and return coil output exhibiting a relatively narrow 
frequency band indicative of that band of the nuclear species. In the wide 
5 band, the contribution of particular spins to the return coil output is not 

dependent upon their location within the specimen. Therefore, the output does 
not indicate the location of components in the specimen. The frequency and 
decay of the output can be used only to identify components of the specimen. 

In order to produce a spatial image of the specimen, gradients are 
10 established in the polarizing field. Although the direction of the polarizing 
field remains the same, its strength varies along the x, y, and z axes oriented 
with respect to the specimen. Varying the strength of the polarizing field 
linearly along the x-axis, causes the Larmor frequency of a particular nuclear 
species to vary linearly as a function of its position along the x-axis. With 
15 magnetic field gradients established along the y-axis and z-axis, the Larmor 
frequency of a particular species will vary linearly as a function of its position 
along these axes. 

By performing a Fourier analysis of the return coil's output, the 
frequency components of the output can be separated. With a narrow band 

20 excitation field applied to excite a select nuclear species, the position of a spin 
relative to the xyz coordinate system can then be determined by assessing the 
difference between the coil output frequency and the Larmor frequency for that 
species. Thus, the MRI system can be constructed to analyze frequency at a 
given point in time to determine the location of spins relative to the magnetic 

25 field gradients and to analyze the decay in frequency to determine the 
composition of the specimen at a particular point. 

Special equential operation of one or more main polarizing field coils, 
polarizing gradient field coils, rf excitation field coils, and return field coils 
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results in generation and sensing of the fields required for proper operation of 
an MRI system. The same coil arrangement can be used to generate the 
excitation field and sense the return field. As described in U.S. Pat. No. 
4,843,322 (Glover); U.S. Pat. No. 4,868,501 (Conolly); and U.S. Pat. No. 
5 4,901 ,020 (Ladebeck et al.), different sequences have been developed for 
specific aspects of MRI system operation. 

Production of angiograms is but one application of conventional MRI 
systems, A variety of pulse sequences and processing techniques have been 
developed for use in MRI angiography [U.S. Pat. No. 4,516,582 (Redington); 

10 U.S. Pat. No. 4,528,985 (Macovski); U.S. Pat. No. 4,647,857 (Taber); U.S. 
Pat. No. 4,714,081 (Dumoulin et al.); U.S. Pat. No. 4,777,957 (Wehrli et al.); 
and U.S. Pat. No. 4,836,209 (Nishimura)]. Blood vessels are readily 
differentiated from surrounding tissue by the pulsatile flow of blood. For 
example, if the excitation field is pulsed at systole and diastole, the contribution 

15 of blood flow to the return field will differ, while the contribution of static 
tissue and bone to the return field will be the same. The contribution from the 
blood vessel is determined by subtracting one return from the other, canceling 
the static component. 

Because neuronal tissue does not exhibit the flow-distinctiveness of 
20 blood vessels, MRI angiography systems and pulse sequences cannot be used 
to generate suitable images and conventional MRI systems and sequences used 
for general imaging of tissue and bone do not provide acceptable results. 

Contrasting agents 

The use of pharmaceutical agents to enhance the contrast of neural tissue 
25 relative to surrounding tissue in the images produced is one technique proposed 
for use in enhancing the imaging of neural tissue [e.g. PCT EP 91/01780 (Filler 
et al., WO 92/04916, 1992)]. Therein, a two-part contrast agent was injected, 
taken up, and transported, by the nerve of interest. The first part of the agent 
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promoted neural uptake; the second the desired image. Injected into muscle, 
the agent undergoes axoplasmic flow in the nerve supplying that muscle, 
thereby tagging the nerve. The second part of the agent has a magnetically 
active component. 

5 Because of the ability to image only a single nerve or nerve group, an 

increasing preference to avoid the use of invasive procedures whenever 
possible, and the typical reduction of the intensity of the imaged nerve, the use 
of contrast agents has certain limitations in neural imaging. 

Anisotropy 

1 0 MRI has been used somewhat successfully without contrast agents to 

map white matter nerve tracts in the brain. White matter coursing through gray 
matter tissue in the brain , unlike the surrounding gray matter, exhibits 
relatively high anisotropic diffusion. Moreover, in axonal pathways 
surrounded by myelin sheaths, water mobility is relatively high, while water 

1 5 mobility perpendicular to the tracts is low. [Douek et al., Myelin Fiber 
Orientation Color Mapping, BOOK OF ABSTRACTS, SOCIETY OF 
MAGNETIC RESONANCE DSf MEDICINE, p. 910 (1991)] 

Briefly, a plurality of perpendicular and parallel field gradient pulses 
(diffusion gradients), oriented relative to the white matter tracts to be imaged is 

20 used. Pulsed gradients change the received signal phase from all of the spins. 
The relative effect of these diffusion gradients is the canceling out of stationary 
spins. Spins moving from one spatial position to another in the time between 
the two diffusion gradients, on the other hand, experience changes in the 
frequency and phase of the spin magnetization. The net effect is a reduction in 

25 the received signal, which is greatest for spins diffusing the farthest between 
the two pulsed gradients. 

Diffusion-weighted imaging 
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Thus, MRI is one of the most versatile non-invasive imaging techniques 
in biomedicine, providing functional as well as anatomical information. MRI 
can also provide images exhibiting quantitative information regarding motion 
of water molecules. Given the anisotropic nature of white matter, water will 

5 diffuse freely along a tract, but is restricted in its motion perpendicular to the 
tract. When the diffusion gradient is aligned with the tract, there is thus a 
greater reduction in signal than when the diffusion gradient is aligned 
perpendicular to the tract. Because this phenomenon is not exhibited by the 
surrounding gray matter tissue, the white matter tracts can be identified. MRI 

10 dealing with microscopic motion within a single voxel is referred to as 
diffusion weighted imaging (DWI). 

The MRI sequence for DWI, where differences in intravoxel water 
motion (apparent diffusion coefficient, D app ) are treated as another contrast 
mechanism. This is an extension of the original Stejskal-Tanner sequence for 
1 5 measurement of molecular diffusion by nuclear magnetic resonance (NMR), 
which employed diffusion gradient pulses (time-dependent field gradient) to 
encode quantitative information for molecular motion (diffusion coefficient) 
into a signal intensity. 

As noted above, intravoxel anisotropic water motion produces signal 
20 attenuation of DWIs dependent on the direction of the diffusion gradient pulses 
applied. This directional dependency, most conspicuous in the myelinated 
fibers, is likely due to anisotropic restriction of water diffusion. In contrast to 
conventional DWIs, where diffusion gradient pulses are applied to three axial 
directions simultaneously, DWIs obtained using diffusion gradient pulses 
25 applied to only one spatial axis are generally referred to as anisotropic. 



While intravoxel isotropic water motion produces identical effects on 
signal intensity of anisotropic DWIs regardless of the direction of the diffusion 
gradient pulses, intravoxel anisotropic water motion affects signal intensity in a 
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manner dependent on the angle of the direction of anisotropy with respect to 
the spatial axis in which the diffusion gradient pulses have been applied. 
Therefore, each anisotropic DWI can be considered the projection image of 
intravoxel anisotropic water motion. Each pixel of anisotropic DWIS contains 
5 the spatial information in three-dimensional resolution, but carries only one 
dimensional information regarding anisotropic water motion in space, resulting 
in image resolution that is inferior to resolution obtained where each pixel has 
three-dimensional information. 

High angular resolution diffusion-weighted MRI 

1 0 The sensitivity of the magnetic resonance (MR) signal to molecular 

diffusion provides the most sensitive non-invasive method for the measurement 
of local tissue diffusion characteristics (1). The basic effect from which 
diffusion information is derived is the signal diminution due to diffusive 
motions along the direction of an applied gradient field (2, 3), The fact that 

1 5 diffusion can have a directional dependence was recognized early on (4) but 
found a great resurgence of interest in the application to diffusion weighted 
imaging of human tissues, where inferences about tissue structure can be made 
from the directional dependence it imposes on the local diffusion. Anisotropic 
diffusion was first demonstrated in the brain by Moseley (5, 6) and has been 

20 used to study a variety of other tissues (7). The determination of anisotropy 
requires a reconstruction of the local apparent diffusion, D app . If the diffusion 
is process has spatially homogeneous Gaussian increments, the directional 
dependence can be completely characterized by the diffusion tensor D that 
relates the signal loss along an applied gradient in an orthogonal Cartesian 

25 system defined by the imaging coordinate system to the diffusion along a 
direction in an arbitrarily rotated orthogonal system defined by the tissue (8). 
Reconstruction of local Gaussian diffusion can thus be posed as one of 
estimating the diffusion tensor (8), which, in principle, requires only 6 
measurements plus an additional measurement for normalization (9). This 
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technique is of particular interest in its application to the characterization of 
white matter tracts (10, 1 1). 

However, it has long been recognized that the Gaussian model for 
diffusion can be inappropriate within the complex structure of human tissue (8, 
5 12). One way in which this model can fail is the presence of multiple fiber 
directions within a single imaging voxel. Because the diffusion tensor model is 
no longer appropriate, the characterization of the diffusion in such voxels 
becomes problematic. A novel approach to this problem, proposed by Tuch 
(13), is to map D app at high angular resolution in order to more accurately 
10 detect variations in diffusion along different directions. This has been extended 
to a scheme by Wedeen (14) where measurements through a range of diffusion 
sensitivities are made. There remains, however, no method for characterizing 
the diffusion measured by these high angular resolution diffusion weighted 
(HARD) methods. 

1 5 In view of the above discussion, it is evident that there exists a strong 

need for an effective method of reconstructing and characterizing three 
dimensional information from axial anisotropic DWI projection images. The 
present invention fulfills these needs and provides additional advantages that 
will be apparent from the detailed description hereinbelow. 

20 

Summary of the Invention 

Accordingly, an object of this invention is to provide a method and 
system for the computerized analysis of MRI data especially as it relates to 
visualizing white matter properties and/or structure. 

25 The invention contemplates a method for identifying diffusion 

anisotropy without invoking the diffusion tensor formalism. The method 
comprises the collecting of a plurality of high angle resolution diffusion image 
data employing a diffusion-weighted stimulated echo spiral acquisition process. 
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The method further contemplates computing the spherical diffusion variance of 
the diffusion data in each voxel by a spherical harmonic transform algorithm, 
identifying components in the three diffusion channels, wherein diffusion 
channels are broken down into direct sum subspaces representing isotropic, 
5 single fiber, and multiple fiber components, and wherein asymmetries 
produced by artifacts fall into channels impossible to reach by diffusion, 
thereby, providing a direct means of noise reduction within said channels as 
well as means for identifying artifactual effects. The method is also useful for 
determining magnitude and direction of diffusion by computing from the 
1 0 transform magnitude and phase. 

Also contemplated is a method for characterizing multi-component MRI 
images of multidirectional crossed fibers in a specimen by obtaining a plurality 
of high angular diffusion weighted images, each image obtained by applying a 
diffusion gradient pulse, and employing a simple spherical harmonic transform 
15 algorithm for identification of diffusion anisotropy based upon the variance of 
the estimated apparent diffusion coefficient as a function of measurement 
direction. The specimen can be white matter in the mammalian body. 

Also contemplated is a spherical harmonic transform algorithm useful in 
characterizing diffusion anisotropy from measurement data collected by 
20 employing high angular resolution magnetic resonance imaging, said transform 
derived from mathematical group theory. 

Further contemplated is computer software programmed to perform 
rapid electronic characterization of MRI images employing group theory and 
the spherical harmonic transform as described herein below. 

Brief Description of the Drawings 

The foregoing summary, as well as the following detailed description of 
the invention, will be better understood when read in conjunction with the 
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appended drawings. For the purpose of illustrating the invention, there are 
shown in the drawings, certain embodiment(s) which are presently preferred. It 
should be understood, however, that the invention is not limited to the precise 
arrangements and instrumentalities shown. 

5 Figure 1 diagrams a typical spherical coordinate system used in physics. 

Figure 2 shows the mean squared error between the true diffusion log (signal) 
(Eqn 13) and the approximate form (Eqn 14) as a function of 6- values and 
azimuthal angle between two identical fibers oriented in the x -y plane (ie, 
polar angle 0 = 90°). The approximation error gets worse for large b and l^ge 
10 but is less than .05 for b < 5000*/^. 

Figure 3 depicts the effect of the coupling term in the approximation Eqn 14 in 
a voxel containing two fibers oriented in the x -y plane at an angle ta<j> = 90° 
relative to one another for three 6-values: b = 500, 2500, 5000 9haaa . (Left) The 
true D app and the approximate D m determined from Eqn 14. (Right) The 
15 mean and coupling terms in Eqn 14. For increasing 6-values the approximation 
is less accurate, as demonstrated by the increasing discrepancy between the true 
and the estimated shapes in Figure 2. 

Figure 4 is a graphical representation of the Spherical Harmonic 
Decomposition. The contributions to D app from the different L orders of the 
20 spherical harmonic decomposition are shown for a single fiber (a) and a double 
fiber (b) voxel. 

(a) Single fiber. The D app contains contributions only from the L = 0 term (the 
sphere) and a more complex shape produced by the L = 2 spherical harmonic 
components. 

25 (b) The D m for two identical fibers oriented * relative to one another, contains 
contributions only from the L = 0 and L = 2 spherical harmonic components, as 
in (a), as well as a contribution from the L = 4 component (bottom right). Note 
that the four-lobe structure of D app in the equatorial plane is generated by the 
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addition of the 4 positive red lobes and the 4 negative blue lobes of the L = 4 
component to the L = 0 component sphere. The color scheme for the 
coefficients has red at maximum positive and blue at maximum negative, while 
for D app is max positive (red) to zero (blue). 

5 Figure 5 is a graphical representation of the spherical harmonic transform. For 
each value of L, there corresponds 2L + 1 values of M:-L, ... ,0,. . . >L. The box 
at each (L,M) coordinate contains an array (upper right diagram) whose 
coordinates are varied parameters at the fixed values (L,M). The transform is 
identically zero for \m\ > L (blank boxes). 

10 Figure 6 shows rotational variations of the SHD. Reconstructed D app s are on 
the left. The SHT coefficients are on the right with the varied parameters 
azimuthal (A<fi) angle on the horizontal axis and the polar (A0) angle on the 
vertical axis and displayed as described in Figure 5. (a) Single fiber rotated 
through the full range of (# 0). (b) Two identical fibers oriented (A# A$) 

1 5 relative to one another, (b left) D app for (A3 A$ = (90°,0°). (b right) Spherical 
harmonic transform of two fibers with one fiber rotated through the full range 
of (A 6! A0). The angular variations are with respect to the first fiber which is 
fixed along x = 0. The SHT in (a) and (b) shows that only L = 0 9 2 components 
arise from the single fiber and only L = 0, 2, 4 arise from the multiple fiber. 

20 This is a consequence of the fact that rotations only mix M components. 

Figure 7 shows that MR diffusion measurements are sensitive only to the 
absolute value of the direction: motion in x has the same effect as motion in -x. 
This imposes a projective symmetry on the measurements so that antipodes on 
the measurement sphere (Top) are identical. (Middle b-d) Projective subspace 
25 of a single fiber. (Bottom e-g) Projective subspace behavior of the coupling 
term. 

(Top) Antipodes on the measurement sphere (red dots). 
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(Middle)(a) Single fiber in the equatorial plane (9 = 90°). Only the M= 0, 2 
components of the L = 2 channel contribute, (b) Single fiber oriented along the 
x axis (9 = 0°). Only the M = 0 component of the 1=2 channel contributes, 
(c) Single fiber oriented at (9 = 45°, ^ = 0°). All of the Mcomponents of the L 
5 =2 channel contribute. (Bottom)(d) SHD components for two identical fibers 
oriented at A9, = 90° in the meridian plane (0 = 90°). Only the even \m\ <2 
components of the even L components contribute, (e) SHD components for 
two identical fibers oriented at A# = 90° in the equatorial plane. Only the M= 
0, 4 components of the L - 4 channel contribute, (f) SHD components of the 
1 0 coupling term for two identical fibers oriented at y = 90° with the first fiber 
oriented along 0) = (45°, 90°) and the second tilted second in the meridian 
plane (<f>, 6) = (135°, 45°). All M components of the even L components now 
contribute. 

Figure 8 depicts single fiber rotations. (Top a, b) Single fiber azimuthal 
1 5 rotation. Rotations in the equatorial plane produce only phase changes. 
(Bottom c, d) Single fiber polar rotation. Rotations in the meridian plane 
produce only magnitude changes. 

(a) Single fiber fixed in the equatorial plane (0= 90°), at three different 
azimuthal angles of rotation: <f> = (0°, 45°, 90°). The magnitudes of the 

20 components remains unchanged, but the phase is altered. 

(b) The phase ^ / 2 of the coefficient a, 2 , 2 of the SHT shown for the fiber in the 
equatorial plane (0 = 7tl 2) rotated through the range of orientations in the 
equatorial plane 0 < <j> < 1 80*\ 

(c) Single fiber fixed in the meridian plane (<f> = 90°), at three different polar 
25 angles of rotation: 0 = (0°, 45°, 90°). The magnitudes of the components are 

redistributed amongst the available L = 2 channels (Af = -2, 0, 2) but the phase 
remains unchanged. (Bright white pixels are numerical noise artifacts aliased 
from 0 to %) 



WO 02/082376 



PCT/US02/12668 



-13- 

(d) The magnitude of the coefficient a 2t 0 of the SHT (solid line) for a fiber in 
the meridian plane = 0) rotated through the range of orientations in the 
meridian plane 0 < 0< 180°. The circles are the function cP 2 (Cos 9) (with a 
scaling factor c) and show exact correspondence with the magnitude variations 
5 oia 2% o. 

Figure 9 is about double fiber rotations. (Top a-b) Double fiber azimuthal 
rotation. Rotations in the equatorial plane produce only phase changes in the 
coupling term. (Bottom c-d) Double fiber polar rotation. Rotations in the 
meridian plane produce only magnitude changes in the coupling term and a 
1 0 - mixing of energy amongst the M components. 

(a) The coupling term for two identical fibers in the equatorial plane (0=%/2) 
at a fixed orientation of 90° relative to one another, rotated relative to the 
laboratory coordinate system through <f> = (0, % 1 1 6, 7t / 8) radians. The shape 
of the coupling term does not change, but its orientation relative to the 

1 5 laboratory changes. As a result, the phase of the coefficients changes, but their 
magnitude stays the same. 

(b) The phase 01 4 of the coupling term for two identical fibers in the 
equatorial plane (0= 7t / 2), at a fixed orientation of 90° relative to one another, 
rotated relative to laboratory coordinate system through a range 0 ^ 0 < % 

20 radians. 

(c) The coupling term for two identical fibers in the meridian plane (0 = 7t / 2) 
at a fixed orientation of 90° relative to one another, rotated relative to the 
laboratory coordinate system through <f> = (0, n 1 16, n 1 8) radians. The shape 
of the coupling term does not change, but its orientation relative to the 

25 laboratory changes. As a result, the relative magnitude of the components with 
L changes, but the phase stays the same. 

(d) The relative magnitude of the SHD components of the coupling term for 
two identical fibers in the meridian plane (^= 0), at a fixed orientation of 90° 
relative to one another, rotated relative to the laboratory coordinate system 
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through a range 0 < 6 <> % 1 2 radians . The angle determined by the changing 
magnitudes in the even and odd M channels is proportional to the orientation: 
<f>/4 = tan 2 (Y 4 °/Y 4 2 ). 

Figure 10 shows double fiber relative rotations and the Fractional Multifiber 
5 Index. 

(a) The coupling term for two identical fibers in the equatorial plane (0= rc / 2) 
at three different relative orientations A<p = (0, n 1 4, k 1 2) radians relative to 
one another, with the mean angle between them fixed relative to the laboratory 
coordinate system $ = (0, tc / 16, % 1 8) radians. The shape of the coupling term 

10 does not change, but its size changes. As a result, the magnitude of coefficients 
of the coupling coefficient grows with increasing relative angle, but the phase 
stays the same. 

(b) The energy in the L = 4 component of the SHD components of the coupling 
term for two identical fibers in the meridian plane (0 = 0), for a range of 

1 5 relative orientations 0 < A < for a fixed mean angle relative to the laboratory 
coordinate system (0°). The maximum magnitude occurs for A = 90°, ie., when 
the fibers are perpendicular. 

(c) (Bottom Left) The Fractional Multiplier Index (FMI) for a simulation of 
two identical fibers for a full range of relative azimuthal A<f> and polar AO 

20 orientations between the fibers, sampled in increments of 7.2° in both A<p and 
AO, (Bottom Right) A few of the corresponding local diffusion D app ((j> t 0). 

Figure 11 shows diffusion encoding directions generated by the pulse sequence 
are spherical tessellations of an icosahedron of user-specified order (shown 
here for clarity is the fifth order tessellation). 

25 Figure 12 demonstrates spherical harmonic transform and the separate fiber 
channels it produces, (a) Compound spherical harmonic transform 
components of HARD data collected on a normal human volunteer with b = 
3000 s/mm 2 and 43 diffusion encoding directions. As predicted, the energy is 



WO 02/082376 



PCTYUS02/12668 



-15- 

indeed confined to the L = 0, 2, 4, except for the eddy current artifacts 
appearing at (L, M) = (0, 1). 

Figure 13 shows the shape of D app estimated from the fiber channels in Figure 
12. (Top Left) Single fiber channel (greyscale) overlayed in color with 
5 multifiber channel (FM > A) (color scale). (Top Right) Greyscale image is 
relative anisotropy index determined from diffusion tensor calculation from 
same data. In circles, there are three points from the isotropic (A), single fiber 
(B) and (C) multiplier channels for which D app is (colored shapes) reconstructed 
from all three SHT components L = 0, 2, 4. The gray matter voxel (A) is 

1 0 essentially isotropic, so that D app is a sphere and unaffected by the inclusion of 
L = 2, 4. The single fiber voxel (B) has the characteristic peanut shape , which 
is unaffected by the inclusion of L = 4. The voxel in (C) requires L = 2, 4 to 
represent D app , which in this case is consistent with two identical but 
perpendicular fibers. The DTI reconstruction (c), on the other hand, can only 

15 accurately reconstruct the isotropic and single fiber channels. Both, (A) 
isotropic and (B) single fiber voxels are evident, but (C) multiplier channel 
appears black, being the extreme case of two fibers crossed at 90°. 

Figure 14 shows fiber magnitude and orientation. (A) Magnitude of single 
fiber channel, representing the presence of single fibers, (B) map of estimated 
20 0, and (C) map of estimated 0. 

Figure 15 shows: (Left) a white matter map reconstructed by adding the energy 
in the L = 2, 4 channels, and (Right) a map made from the energy in the odd 
order components (L = 1, 3, 5) showing little structure. These images are not 
on the same scale, but scaled independently for display. 

25 Figure 16 represents evidence of consistency with the two-fiber model. D app 
measured from multifiber channel in Figure 12 (top) and simulated (bottom) 
for (1-r) // =f 2 ,dt = 90°,f = 2/ 2) d<fi = 90°,f.= 2f 2 , dip = 75°. The shapes have 
been interpolated from the diffusion directional sampling for clarity. 
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Description of the Preferred Embodiment 

MR Diffusion Measurements 

5 Before presenting our strategy for characterizing diffusion anisotropy in 

multifiber systems, we summarize in this section the basic mathematical 
underpinnings of MR diffusion measurements. In what follows, a fiber will be 
defined as a particular diffusion tensor D. With this definition, a large bundle 
of parallel fibers would be synonymous with a single "fiber". While it is 
1 0 natural to confer cylindrical symmetry on a diffusion tensor as part of the 

definition of a fiber, we relax this restriction in order to more clearly establish 
where such symmetry actually effects the more general results. 

Single Fiber 

We first consider diffusion measurements made in a voxel containing a 
15 single fiber, following Hsu and Mori [17] throughout. The signal attenuation 
can be written 

S(b,D)=S 0 e- bDapp 

where b is related to the fc-space trajectory k (t) byb = Jo™ I^(t)dt and 
20 incorporates the gradient strengths. The gradient directions can be defined by 
the unit vector u. If the measurements are made along the principal axes, i.e., 
in the coordinate system in which the diffusion tensor is diagonal, the apparent 
diffusion coefficient D app is ([17]) 

D apP = u T • D A ■ u 

25 

and the applied diffusion encoding gradient u is in the direction of the 
eigenvectors of 2). 
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u= u 2 , i^A= 0 d 2 0 (3) 
V u 3/ \0 0 d 3 j 

where the eigenvalues are the principal diffusivities {d^ dz, d^}. Generally, 
however, this principal axis coordinate system is not known. The applied 
5 diffusion encoding gradients v are therefore not coincident with the principal 
axis system, but are related to it by a rotation R: 

u = Rv (4) 

where v is a unit vector in direction of diffusion encoding. 

Thus, one usually wants to infer the principal diffusivities and the 
10 rotation R. From these can be determined the diffusion properties, such as the 
anisotropy, and the fiber directions. The rotation R is defined within the 
coordinate system shown in Figure 1 . The two angles that define the direction 
in this coordinate system are the polar angle Bz [0, n ] which is defined as 
the angle between the vector and the positive z-axis, and the azimuthal angle <f> 
15 e [0,27t ) , which is defined as the angle in the x - y plane relative to the 

positive x axis. It is also common to use the elevation angle 8 = 9O°-0, which 
is the angle between the vector and the x - y plane . This is often denoted by fa 
however, as in Hsu and Mori [17]. We will retain the standard physics usage, 
depicted in Figure 1, where (0, <fi) denote the polar and azimuthal angles, 
20 respectively. We will often employ the shorthand notation Q s (fi 9 

It is helpful to note that {0, <f> ) are two of the Euler angles [18], used to 
describe rotations in 3-dimensional coordinates. These angles typically 
denoted (a,p,y ) where a is the azimuthal rotation angle, P is the polar rotation 
angle, and y is a rotation about the new axis defined by the rotation through 
25 ( a,P). For the description of a single point (i.e., a measurement) on a sphere, 
as is the case in this paper, rotations about the final (radial) axis are 
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unimportant, so the rotations can be described by the two angles ( a,P ). It is 
common in this case to denote these (6^). 

The gradient direction vectors in the two coordinate systems are related by a 
5 rotation ([17]). 



20 



R = 



f sin# — cos<£ 0 ^ 
cos <f> cos 6 sin <f> cos 9 - sin 6 
^cos^sinfl sin <j> sin 9 cos 6 j 



(5) 



10 The apparent diffusion coefficient for an arbitrary gradient direction vcan thus 
be written ([17]) 

D app = v T D-v (6) 

where 

15 D=zR T D A R (7) 

Eqn 7 defines the diffusion tensor D in a rotated coordinate system. The signal 
from a single fiber is typically expressed in the form 



log(S/S 0 ) = -&£> (8) 

where we introduce the shorthand notation D = D app =v T Dv. For any 
symmetric matrix D, such as the diffusion tensor, the product x T Dx is a pure 
quadratic form [19]. From Eqn 8 one can estimate the diffusion tensor D by an 

25 eigenvalue decomposition whose eigenvectors effectively determine the 
rotation of the fiber coordinate system relative to the laboratory system and 
whose eigenvalues determine the diffusivities. Since/) is positive definite, it 
can be written in the form of Eqn 7 where D A is diagonal and the unit 
eigenvectors of D are the columns of R. The rotation y = Rv produces the sum 

30 of squares 
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n 

v ?Dv = v T RD A R T y = y T D A y = £ A, (9) 



From [19], the equation of \?Dv - 1 describes an ellipsoid whose axes end at 
5 the points where A^ y ? =1 and where the remaining y components are zero. 
Undoing the rotation, these points are in the directions of the eigenvectors and 
the axes have half-length 1/ It is important to emphasize, however, that the 
ellipsoid that describes the eigenspace of the diffusion tensor is not a 
description of the shape of the measured local diffusion [15]. 

10 Multiple Fibers 

The definition of a fiber above can be extended to multiple fibers by 
defining ihek'th fiber as synonymous with the k'th diffusion tensor D k with 
no assumption of cylindrical symmetry until it is necessary. In addition, the 
assumption will be made throughout that there is no exchange between fibers 
15 so that the signals add independently. 

In general, the signal from a voxel containing n fibers can be written 

S = SoX>e-^* (10) 

20 Where / k is the volume fraction of the k'th fiber (iT^ /pi). It will prove 
useful to express this in the following form (Taylor expanding about b = 0): 

log(S/S 0 )« -6 y|EA(l-/ fc )A£)| 1 -2 E /i/;AAiAD,i| (11) 

25 where we have defined the differential apparent diffusion AD/cj = D k -Dj. 
This form of writing the signal is useful in that it expresses the effect of 
additional fibers I&2 relative to the first (k = 1). For a single fiber, the terms 
in the bracket drop out and the log signal assumes the pure quadratic form 
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utilized in Eqn 8 in the standard diffusion tensor approach. Each of the D is a 
quadratic, so each term in Eqn 1 1 is an even power polynomial in gradient 
direction v. The approximation to orders O (b 2 ) yields even polynomials up to 
order 4. 

5 Special case: Two fibers 

For multiple fibers within a voxel, Eqn 1 1 can be complicated. 
However, the special case of two fibers (n = 2) is instructive and a useful 
approximation in practice. This is the simplest multifiber case, and the 
resulting equations have an intuitively clear interpretation and are numerically 
1 0 relatively simple to manage. From Eqn 1 1 , the log signal case be written 

log (S/So) « -6 [f x D x + f 2 D 2 ] + fx h b 2 A£)|i (12) 

where / 1 and/i are the volume fractions in compartments 1 and 2, respectively 
so that f\ +f 2 = 1 and AD 2 \ =Z D 2 -Du and we have kept terms up to second 
order in b. The measurements D i are composed of second order polynomials, 
15 so the coupling term £ s f x f 2 b 2 A D 2 2 i is composed of even order 
polynomials up to fourth order. 

The two terms linear in b represent the individual fiber components and 
are pure quadratic forms. In addition, there is a coupling term, second order in 
b y with coefficient f\ f 2 that is of fourth order. Note the interesting fact that the 

20 coefficient term f x f 2 is a quadratic function, with a maximum at/i - f 2 . The 
magnitude of the coupling thus depends on the volume fractions. However, the 
shape does not. This is clear from the fact that the volume fraction enters only 
as a multiplicative factor in £ Variations in the shape of the coupling term are 
more easily understood by rewriting in a more illuminating form. First note, 

25 from Eqn 6 and 7, that A f) n = v T R\b l2 R l v (13; 



Where we have defined Z>i 2 s RT 2 £>2,\ R u - -°i,a 

R\1 = Jf?2-Ri * 



(14) 
(15) 
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The term R i2 is the product of rotation matrices and is therefore also a rotation 
matrix (by virtue of the group properties of the 3-Dimensional Rotation Group, 
denoted O (3) [20]). It first undoes the rotation R j then applies the rotation R 2 . 
If the rotations R\ and R 2 are the same, meaning that the fibers are pointing in 
5 the same direction, then R x2 = 1 , and D l2 is a diagonal matrix with eigenvalues 
equal to the difference in principal diffusivities between Dy and D 2 . If the 
fibers are identical in orientation and diffusivities (by our definition, D x -D 2 \ 
then D l2 = 0, so that AZ) 12 = 0 and the coupling term disappears, reducing the 
log signal to the standard single fiber form (Eqn 8). In general, though, the 
1 0 term A D 2i is of the same form as the rotated single fiber diffusion matrices 
(e.g., Eqn 7) but in terms of the new tensor D l2 , which we shall term the 
reduced diffusion tensor. 

The accuracy of the approximation in Eqn 12 depends upon both the b- 
value and the relative orientation of the fibers. As shown in Figure 2, this 
15 approximation is good up to quite high 6-values. The manifestation of these 
errors is most clearly visualized by plotting Z> app as a function of 6-value for 
AO- 90° (the angle at which the influence of the coupling term is the greatest), 
an example of which is shown in Figure 3. 

The results for HARD diffusion measurements can be summarized as follows: 
20 1 . The measured local diffusion is related to the diffusion in the fiber 
coordinate system by three-dimensional rotations. 

2. The measured log signal from a single fiber is a pure quadratic form. 

3. The measured log signal from a multiple fiber can be approximated by 
the sum of two pure quadratic forms and a coupling term of even 

25 polynomial up to order 4. 



The question is this: If an eigenvector decomposition is sufficient to 
determine the diffusion tensor in the single fiber case, is there a decomposition 
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sufficient to determine the diffusion in the case of multiple fibers? We show in 
the next section that this is indeed possible. 

The Spherical Harmonic Decomposition 

Spherical Tensor Representation of Rotations 

As described in Section 2, the diffusion measurements along different 
encoding directions can be expressed as rotations in 3-dimensions relative to 
the (unknown) principal axis system of the diffusion. In the HARD technique, 
these measurements are along a set of directions covering the range of the 
spherical coordinates (6J $\ The rotations R were expressed using a form or 
representation in terms of the Euler angles. However, the problem was 
formulated in the familiar Cartesian coordinate system, or basis e = {*, y, z}, 
because this is the natural coordinate system for imaging. Systems with 
spherical symmetry are often more conveniently handled in the spherical basis: 
e -{r, 0> <f>} as defined in Figure 1. As we show next, rotation matrices 
transformed into the spherical basis are the spherical harmonics [21]. Tensors 
transformed into this representation are called spherical tensors [21]. 

The HARD diffusion measurements have an inherent spherical 
symmetry because they are made by a series of three-dimensional rotations. 
The inherent symmetry in this problem can be elucidated through the theory of 
groups, which was originally developed for the purpose of characterizing 
symmetry [21]. The concept and consequences of groups, although extremely 
powerful, are conceptually relatively simple. Their great power is that they 
facilitate the characterization and classification of mathematical structures into 
classes or groups with similar symmetry properties. All members within a 
particular group can then be treated as equivalent, even if their specific 
manifestations differ. For example, two points constrained to move on the 
surface of the same sphere can be seen as having identical symmetry properties 
even if their precise locations on that sphere differ. 
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The importance of the spherical tensor formulation is encapsulated in 
the following expression, which shows how a spherical tensor J is affected by 
a general rotation & m » m {M) in some representation 5R of the rotation group in 

the basis e: 
5 < 

m'=-l 

Notice that rotations of the individual m components within a particular 
component rank 1 do not " mix" elements ainongst components of different 
rank I. This quality of the Tis called irreducibility. Eqn 16 can be considered 
10 the defining quality of a spherical tensor: a tensor that transforms accordingly 
is by definition a spherical tensor. The spherical tensor representation is useful 
because rotations of the individual tensor components preserve the rank. 

Both the spherical tensor and the rotations can be related to our specific 
problem as follows. The diffusion tensor is a second rank tensor and so 
15 consists of nine components, represented by a 3 x 3 matrix, In its irreducible 
representation, the tensor is written as the sum of three terms: 

T = T u + T° + T 5 



20 where 1° = 1° I in which T° is a rank 0 tensor (i.e., a scalar) and / is the 3 x 3 
identity matrix, T* is an anti-symmetric rank 1 tensor *i.e., a vector) and r * is a 
symmetric, rank 2 tensor. The rotations ® l m 'm (a,3,y ) in the spherical tensor 
basis expressed in terms of the Euler angles are called the Wigner rotation 
matrices [21]. For a point in spherical coordinates, the Wigner rotation 

25 matrices are proportional to the spherical harmonics [22]: 
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Therefore, the process of rotating a diffusion tensor can be reformulated by 
expressing the diffusion tensor in an irreducible form in which its individual 
components transform separately under rotations affected by spherical 
harmonic components. In this general formulation, the concept of the diffusion 
tensor can easily be extended to more complex structures by considering 
tensors of higher rank because now the transformation under rotations are of 
exactly the same form (Eqn 16) and the functions that perform the rotations are 
exactly the same basis set (the spherical harmonics). 

A pplication to HARD Measurements 

Now, consider the general case of a HARD measurement of a voxel of 
unknown fiber composition. The measured apparent diffusion coefficient D ( 
Q ) is then an arbitrary real function. The complex spherical harmonics form a 
complete orthonormal basis [18] so an arbitrary real function parameterized by 
the spherical coordinates ( 0, ) can be expanded in a Laplace series: 

^(M)=EEa n (M) 

1=0 m=-l 

The coefficients a lm are determined by [1 8] multiplying both sides of Eqn 19 by 
Y™ '*( 0,<j> )and using the orthogonality condition 

/•27T PIT 

/ / Yr{0,<l>)Y£{6^)sm{0)d0d<l> = & lk 8 mn 
Jo Jo 

The expansion coefficients are uniquely determined by multiplying each side of 
Eqn 19 by Y k n *(0,$) and integrating over the sphere. The result is that the 
coefficient can be determined by 

a im = F* rD app {0, <j>) <t>) sin(0) d9d<j> 
Jo Jo 

This is precisely analogous to a Fourier decomposition of sinusoidal functions, 
but on the unit sphere. This will be called the spherical harmonic transform or 
SHT of the measured apparent diffusion coefficient. 



WO 02/082376 



PCT/US02/12668 



-25 - 

The utility of describing the measured HARD (log) signal in terms of 
group theoretical constructs can now be shown directly as follows. The 
measured HARD signal D (Q) is an arbitrary real function and so can be 
5 expanded in term of a Laplace series (Eqn 1 9) with the coefficients determined 
by the SHD (Eqn 21). But the symmetry of HARD diffusion measurements 
imposes severe restrictions on the expansion coefficients that allow direct 
classification of the diffusion from the SHD. In particular, the SHD of Z)(Q) 
produces the coefficients of the irreducible representation of D(Q) . 

10 Specifically, the following is true: Isotropic diffusion is independent of 

direction, so the lowest order Y 0 o(Q ) is the basis for the representation of D 0 . 
This fiis easily seen because Y 0 o(£l ) is just a sphere, so the calculated 
coefficient merely scales the radius of the sphere. This provides a 2L + 1 = 1 
dimensional representation of the 0 (3). For a single fiber, the irreducible 

1 5 representation of D (Q) (Eqn 1 7) provides a basis for a six-dimensional 

representation of O (3) (with 7* providing 2*0+1 = 1 and V providing 2*2+1 = 
5). The basis functions are then the spherical harmonics Y im (Q) of order L = 
0,2. Because this is an irreducible representation, Eqn 16 expresses the fact 
that fiber rotations do not alter the basis functions. That is, they only produce a 

20 redistribution of energy amongst the M components or variations in the phase 
of the M components, but not the order L. For multiple fibers, the log signal 
from multiple fibers can be expressed in terms of an expansion in even power 
polynomials (Eqn 1 1) and approximated up to relatively high i-values by 
keeping terms up to second order in b, which corresponds to terms of 4' th 

25 order polynomials. A multifiber voxel therefore provides a basis for a 2* (0 + 
2 + 4) + 3 = 15 dimensional representation of 0(3) with basis functions being 
the spherical harmonics of even orders up to order L = 4. The dimensions of 
the representation are the number of measurements that are required to 
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characterize the D (Q). It is for this reason that 2 * (0 + 2) + 2 = 6 
measurements are required to characterize the standard single fiber diffusion 
tensor. 

5 The results can be summarized by the following rather remarkable 

conclusions: 

1 . Isotropic diffusion is described by A 0£? Y 0 °(Q) 

2. Single fiber diffusion is described by Ei^oja lm Y t m {C£). 

3. Multiple fiber diffusion is approximately described by 

10 2/ = 0,2,4 2 maim?™ (O) 

4. In general, the local diffusion, including the magnitude and orientation, 
can be described by sum of spherical harmonics of even order, i.e. 
Ji^atjrr^ / = 0,2,4,... 

5. Odd orders of the spherical harmonics describe asymmetric components 
1 5 and therefore imaging artifacts. 

6. The coefficient aim is determined by the spherical harmonic transform 
of D (Q) 

7. The dimension of the representation is the number of measurements 
required to characterize the apparent diffusion coefficient D(Ci) . 

20 It is important to point out that the order L required to characterize the 

diffusion in a multifiber voxel depends upon the orientation of the two fibers. 
For fibers more closely aligned, higher orders of L will be required. This can 
be seen by considering the simplest case of two identical fibers lying in the 0 
= 90° plane, oriented nearly parallel to one another. Distinguishing between 

25 the two fibers requires high resolution in the azimuthal angle <f> . Since the 
azimuthal dependence of the spherical harmonics is proportional to exp( im<f> ), 
higher frequency variations in tj> require larger values of m and thus higher 
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orders L in the basis functions. We have focused on the simplest case in the 
present paper where L = 4 is considered sufficient, but this is not required. 

The group theory arguments provide a simple and concise description of 
the categorization of voxel diffusion characteristics, since the above results are 
5 expressible as direct sum subspaces. Let us call the " state " of a voxel with k 
fibers y k, where k = 0 means isotropic diffusion. Then we can write 

1 . Isotropic diffusion: V>o = ® i0) - 

2. Single fiber diffusion: * = £ (0) © 9 (2) . 

3. Multiple fiber diffusion: ^2«2> (0) ®2)( 2 >©3)( 4 >. 

10 where the <£(;) are the irreducible representations of the rotation group. In 
principle, the composition of a voxel in terms of these can be determined, as 
well as the magnitude and orientation of the local diffusion. 

This theory is easily confirmed by simulation, as shown in Figure 4 
where the D app for a single fiber and for two identical perpendicular fibers are 

1 5 broken down in terms of the separate contributions from the different spherical 
harmonic components. In these simulation we use a fiber model that possesses 
cylindrical symmetry, since this is a reasonable physical model for white matter 
geometry (e.g., [17]) though the above results do not require this. The 
corresponding SHD are shown in Figure 6, which shows on the right the SHT 

20 of the D app from a single fiber (top) rotated through the entire range of both 9 
and ^, and two fibers (bottom) with one fixed along the x-axis and the relative 
angle between the fibers rotated through the entire range of both 0and <j>. The 
D Qpp for .one particular orientation is shown on the left. No energy is produced 
in channels other than those predicted by the group arguments above. 

25 The symmetry inherent in this problem precludes the power in the odd L 

channels. Energy in these channels in the SHD of actual experimental data is 
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therefore produced by non-diffusion effects, such as slice offsets or eddy 
currents, which are not constrained to the same symmetry properties. This can 
be used to advantage as a means of identifying non-diffusion related variations. 
Of course, such variations may also have power in the diffusion related 
5 channels, so their reconstruction may not be trivial. The information in the odd 
channel might also be useful in incorporating eddy current correction into the 
image reconstruction. 

The ability to characterize the diffusion does not imply that extraction of 
fiber information is easy, however. This becomes apparent in examining the 
10 coupling term £ , about which can be said: 

1 . The shape of £ depends on the eigenvalues of the reduced 
diffusion 

tensor, and thus on the relative anisotropics of two fibers. 

2. The orientation (and hence phase) of £, depends upon the mean 
1 5 orientation of fibers. 

3 . The magnitude of £ depends on the volume fractions and the 
relative orientations of the fibers. 

The last item underscores a basic ambiguity in the diffusion 
measurements: the volume fractions and relative orientations can confound 
20 information about one another in the measured signal. However, the 
orientation affects the phase, whereas volume fraction changes do not. 

Computation and display of the SHD 

Numerical methods. The spherical harmonic decomposition is achieved 
by computing the spherical harmonic transform of the measured (i.e., apparent) 
25 diffusion coefficient (Eqn 21). Unfortunately, unlike the Discrete Fourier 

Transform (DFT), for which there exists a matrix decomposition that allows the 
Fast Fourier Transform (FFT), no such algorithm exists for the SHT. Although 
a variety of algorithms have been proposed for the computation of spherical 
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harmonic coefficients ([23, 24]), no clear algorithm has emerged as clearly 
superior, and the subject remains an area of active research. Therefore, for the 
present work the coefficients were determined by direct computation of Eqn 

21. 

5 The direct computation of Eqn 2 1 on a grid of N values of <f> and M 

values of 0 requires the computation of I x (2L + 1) spherical harmonics 
evaluated at N Appoints. The spherical harmonics involve multiple expensive 
trigonometric evaluations. Many of the trigonometric evaluations are 
redundant, however, so precomputation of these values can be used to speed up 

10 the computation. Direct computation of the measure sin (0) d6d<f> requires M 
trigonometric evaluations and NAfi multiplies. However, an efficient algorithm 
for computing the measures was developed by noting that these weights in the 
summation that approximates the integral are equal to the Voronoi areas for the 
sampling points on the unit sphere. Precomputation of the trigonometric 

1 5 functions and the weights therefore allowed an efficient SHT. This algorithm 
has been incorporated into the author s diffusion plugin module in AFNI [25] 
and is sufficiently fast to rapidly process reasonably large HARD data sets 2 . 
For example, the data shown in the present paper are comprised of 10 slices 
and 43 directions. The SHT on the entire data set took only 1 .25 minutes on an 

20 SGI Octane2 with dual Rl 2000 processors. 

Display of the SHD 

Before proceeding, it is important to outline our basic method of 
displaying the SHD data. Because there is a great deal of information produced 
by this method, we have found it is important to have a concise method for 
25 displaying the results. The SHD calculates components up to a user specified 
value of 4 for all the 21 + 1 values of M associated with each L (because M = - 
L 9 ...0, ... + L). A useful way to display the coefficients is thus on an array of 
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coordinates (L, M) 9 as shown in Figure 5. Since M < L, the diagram has a 
characteristic triangular shape (white boxes have no coefficients). 

Generally, the coefficients are determined for a range of parameters, in 
which case the boxes become arrays wherein the parameters are varied. The 
5 first is where the parameters are the spatial coordinates (x, y). Then each " 
box" on the (L, M) plot is just an image of the spatial distribution of the 
amplitude of that specific (L, M) component of the SHD (e.g. Figure 12 below). 
The other case of importance is when the parameters are the azimuthal and 
polar fiber angles {<f>,0) as shown in the example in Figure 6. This is useful, 
10 for example, in showing that as a single fiber is rotated arbitrarily, there is a 
redistribution of amplitudes within the L=2 channel, but the energy remains 
completely contained with this channel, as predicted. 

The Structure of the SHD of Diffusion Measurements 

Determining significance of multiple fiber channel 

15 The categorization of voxel fiber composition outlined above suggests a 

strategy for the analysis of high angular resolution diffusion data. The 
spherical harmonic transform is taken for each voxel in the image and sorted 
into even and odd order L. The odd orders represent artifacts and therefore can 
be eliminated from the analysis. The remaining even orders up to order L = 4 

20 are then sorted in the following manner. Voxels with significant power in Z=4 
are classified as " multiple fibers", voxels with significant power in L = 2 but 
not L = 4 are classified as single fibers , and voxels with significant power in L 
= 0 but not L = 2 or L = 4 are classified as 44 isotropic" . Although the method 
to determine what is significant involves some subtlety, the basic strategy is 

25 clear. Having determined into which classification each voxel falls, the voxel 
local diffusion is determined from the sum of the appropriate spherical 
harmonics components (isotropic: 1 = 0, single fiber: L = 0,2, multiple fiber: L 



WO 02/082376 



PCT/US02/12668 



10 



-31 - 

= 0,2,4), with the coefficients determined from the spherical harmonic 
transform. 

One method for the determination of significance is suggested by the 
results above that show that the magnitude of the L = 4 term increases with 
increasing relative orientation of the fibers. Therefore a comparison of energy 
in the L = 4 to the L = 2 channel can be used to gauge of whether or not a voxel 
is of state tp 2 or One measure of the significance of a multiple fiber 
channel is the fractional even order greater than 0 in that channel We can 
define the Fractional Multifiber Index or FMI as 

, Leven (22) 



An example is shown in Figure 10 (Bottom) for the simplest example of two 
fibers at a range of relative orientations ranging from parallel to perpendicular. 
This is a reasonable measure of comparison and means of separating single and 
1 5 multifiber channels, begging the question of which threshold to choose. We 
will use this here in lieu of a more complete probabilistic model that will be 
pursued in the future. 

The Symmetry of Diffusion Weighting: Projective Subspaces 

Decomposition of the diffusion into separate isotropic, single, and 
20 multifiber channels is a consequence of the group algebra, which generates 
additive subspaces that depends upon the degree of polynomial necessary to 
describe the measurements. However, this is far from the complete story, for it 
is just the L story. There is structure within the SHD contained in the way in 
which the energy is distributed amongst the M components. 

25 There is a fundamental symmetry in the D app imposed by the imaging 

process because the signal loss due to diffusion along the direction of a 
gradient is insensitive to the sign of the motion. That is, equal diffusive motion 
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in both the +x and -x directions produce the same diffusion related signal loss. 
This results in a projective symmetry, which can be visualized by the diagram 
in Figure 7(top:a). In the mathematics literature, projective is synonymous 
with antipodal, and is unrelated with projection in the sense commonly used in 
5 the physics literature (i.e., the component along a chosen axis). Diffusion 
measurements are therefore represented by the projective subgroup of 0(3) 
(denoted PSO(3) = 0(3)/Z(2) where " / " can be thought of as " mod " and Z 
(2) represents the group of two integers). This means that two antipodal points 
(the two integers) on a sphere are indistinguishable. The effect of projective 
10 symmetry is that it restricts the energy to even values of M because odd values 
are not symmetric in <j> . However, this symmetry imposed on the PS0(2) 
subgroup of 0(3) is only evident for a single fiber in the equatorial plane (0= 
90° ) because arbitrary rotations can be represented by mixtures of (<f>,Q) 
components. These effects are demonstrated by simulation in Figure 7(b-d). 

15 The SHD for two fibers exhibits a similar symmetry. For two fibers in 

the equatorial plane oriented at A0 = 90° to one another, the mean component 
is cylindrically symmetric and therefore does not possess (L, M) = (2, ± 2) 
components. For two fibers in the meridian plane oriented at A 0=90° to one 
another, the coupling term generates all even M terms for L = 0,2,4 except of 

20 (L, M) = (L, ±4). Two fibers in the equatorial plane generates only one set of 
non-zero M components: (L, M) = (4, ±4) components, similar to the single 
fiber, and by symmetry do not possess (L, M) = (2, ±2) components. Again, 
arbitrary rotations produce mixing into all available M components. The 
projective subspace behavior of two identical crossed fibers is shown in Figure 

25 7(e-g). 

Fiber Orientation 



Characterizing the local diffusion amounts to determining difiiision 
tensor(s) and their orientation(s) relative to the laboratory system. In the 
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special case of cylindrically symmetric diffusion, a natural definition of the 
fiber orientation is the direction of the principal (i.e., largest) eigenvector e x of 
the diffusion tensor [8]. For the problem of multiple fiber voxels, the issue of 
orientation becomes more complicated, for one can ask not only the orientation 
5 of the individual fibers, but their orientation relative to one another. In this 
section we reconsider the case of the single fiber, but now in the context of the 
SHD. We then consider the simplest multifiber case of two cylindrically 
symmetric fibers at arbitrary orientations. 

Single fiber orientation from the SHD. 

10 Because the SHD results for a single fiber identically reproduce those of 

the diffusion tensor, a single fiber orientation can always be determined from 
the SHD by transforming from the spherical basis back to the Cartesian basis, 
thereby reconstructing the diffusion tensor, and then determining the fiber 
orientation from the principal eigenvector. 

1 5 The key to determining the fiber orientation from the SHD is to 

recognize that a rotation R<f> of the fiber by an azimuthal angle ^modulates that 
phase of the SHT components because the <f> dependence of the coefficients is 
of the form exp(i ^ ). On the other hand, fiber rotations R e by a polar angle 
produce amplitude variations because the 9 dependence is proportional to Pi 

20 (cos 0), the Legendre polynomial, which is a polynomial in cos 0 of order /. As 
a consequence of Eqn 16, the 0 variations mix energy amongst the available M 
components for a particular L component, but do not exchange energy amongst 
the L components. 

These variations are illustrated in Figure 8(a,b), where a single 
25 cylindrically symmetric anisotropic fiber is rotated through azimuthal angles 
<f> = { 0°,45°, 90° } while fixed at the equatorial plane (i.e. 0= 90° ) and 
through polar angles 0= { 0°, 45°, 90°} while fixed at the prime meridian (i.e. 
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$ = 0 °). The phase and magnitude for the single fiber case shown in Figure 
8(a,b) and Figure 8(c,d) demonstrate that the orientation can be determined 
from the phase, while the 0 orientation can be determined from the amplitude. 

In this example, the phase (Figure 8(a,b)) is determined from the 
5 coefficient a 2> 2 that corresponds to the spherical harmonic Y 2 2 = Vl5/327C 

sin 2 (0) e' 2 * so that the estimate of the azimuthal orientation is ^ = arg [a 2f 2 ]/2 
where arg denotes the phase angle. The magnitude variation of a 2 , o as a 
function of 0 < 0< n is seen (Figure 8(c,d)) to follow P\(cos6) so that the 
estimate of the polar orientation 0 can be made by noting that the a 2> i 
10 component has no energy for 0 =90°, so that the ratio of a 2i i to a 2 ,o can be 
related to the polar angle by: 6 = tan" 1 ( a 2 ya 2 ,i). 

Two-fiber orientation from the SHD 

Determination of the fiber orientation in a multifiber voxel is 
complicated by the fact that the measured signal is not just a function of the 
15 individual fiber orientations, but of their relative orientations and their volume 
fractions. However, the coupling term £ that generates the L = 4 components 
have some remarkable properties that make this problem tractable, at least for 
the two fiber case. Three limiting cases illustrate these properties. 

The first is two fibers coplanar in the equatorial plane (i.e., 0 =90°). If 
20 the angle between the two fibers is kept constant but the two fibers are rotated 
relative to the laboratory coordinate system, then the coupling term £ rotates, 
causing a phase change in the L = 4 components, but no magnitude change. 
This is shown in Figure 9(a-b). 

The second is two fibers coplanar in the meridian plane (i.e., # =0°). If 
25 the angle between the two fibers is kept constant but rotate the two fibers 
relative to the laboratory coordinate system, then the coupling term £ rotates, 
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causing a phase change in the L = 4 components, but no magnitude change. 
This is shown in Figure 9(c-d). 

The third is two fibers coplanar in the meridian plane (i.e., <j> = 0°). If 
the angle between the two fibers is varied but the mean angle between the 
5 fibers is kept fixed, the size of the coupling term £ changes, disappearing when 
the fibers are aligned since this is identical to a single fiber. This is shown in 
Figure 10(top:a-b). 

The orientation results for the special case of two fibers can be 
summarized as follows: 
1 0 • Azimuthal rotations R # of two fibers together (i.e., fixed relative 

orientation but variable mean <f> orientation relative to laboratory frame) 
produces only phase change in the components. 

• Polar rotations R e of two fibers together (i.e., fixed relative orientation 
but variable mean 0 orientation relative to laboratory frame) produce a 

1 5 redistribution of magnitude change in components, but no phase 

changes. 

♦ Variable relative orientation but fixed mean, orientation relative to the 
laboratory frame produces only magnitude changes in the coupling 
component, but the relative amplitudes of the components (i.e. the 

20 pattern) remains unchanged. Utilizing these results necessitates the 

assumption of a two-fiber model. 

Methods 

Images were acquired on a GE SIGNA 1.5T Clinical Imager with high- 
speed gradient hardware using a stimulated echo sequence with a stimulated 
25 echo spiral acquisition previously described [15]. Diffusion sensitive images 
were acquired on five normal human subjects, with approval from the Humans 
Subject Committee at UC San Diego. 
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High angular resolution diffusion encoding is achieved by generating 
gradient directions equally spaced on a sphere by tessellations of an 
icosahedron [13, 15], as shown in Figure 1 1 . This procedure produces 
directions that are equally separated in angle on the surface of a sphere. Single 
5 shot images were acquired at 9 slices with the following parameters: FOV= 
24cm, slice thickness = 3.8 mm, and matrix size 64 x 64 for approximately 
(3.75 mm x3.75 mm x 3.8 mm) isotropic resolution, TR = 2700 ms, TE = 52 
ms. The difiiision parameters were: diffusion gradient duration, 8 = 20 ms, 
stimulated echo mixing time TM- 200 ms, and b = 3000 s/mm , and 43 
10 diffusion directions determined by the icosahedral tessellations of a sphere. 
Twenty averages at each diffusion direction were collected to ensure high 
signal-to-noise ratios, and resulted in a total scan time of = 34 min. 

Results 

The spherical harmonic decomposition (up to order L = 5) of a HARD 
1 5 data set collected in a normal human volunteer is shown Figure 1 2. Note that 
the energy is indeed confined to the L = 0,2,4, except for the artifacts appearing 
at (Z,, M) = (0, 1). This is discussed below. The isotropic, single fiber, and 
multiple fiber channels reconstructed from this transform are shown in Figure 
12 (bottom). Measured D app (0,# ) from the multifiber channel of Figure 12 
20 are compared with simulations in Figure 1 6 assuming 2 fibers oriented 

relative to one another. As predicted, the isotropic channel appears to consist 
of gray matter, the single fiber channel looks like a white matter map, and the 
multifiber channel corresponds anatomically to regions of complex fiber 
geometry. From the coefficients of the SHT shown in Figure 12, the local D app 
25 can be calculated. In Figure 13 are shown representative shapes from the three 
different channels. 

For the representative voxel from the multiple fiber voxel the extreme 
case of nearly identical fibers crossed at 90° was chosen since the failure of the 
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DTI method is most apparent in this circumstance. However, it is important to 
note that regions of significant energy in the multifiber channel do not 
necessarily correspond to " black holes" in the standard DTI maps (such as 
those shown in [15]) since this effect is produced only in the specific case of 
5 identical fibers oriented at 90° to one another. However, these regions do 
correspond to regions in which characterization by a single diffusion tensor is 
incorrect, and will produce spurious results in the estimation of diffiisivities 
and orientations. 

Figure 12 (bottom) shows significant energy in the (L, M) = (0, 1) 
10 channel. Since this is an odd channel, it cannot be related to diffusion. The 
structure of this channel can be understood from the reconstruction of the shape 
produced by the (L, M) = (0, 1). The addition of Y i° to Y 0 ° causes the isotropic 
sphere to be offset in the z-direction. This is consistent with both eddy currents 
and slight imperfections in slice refocusing. The likeness of the Y 0 l channel to 
1 5 a gray matter map is a consequence of the fact that the artifacts that generate 
these components are most evident in regions of large isotropic diffusion where 
the diffusion "sphere" is prominent. 

A white matter map is reconstructed from the energy in the L = 2, 4 
channels, regardless of any categorization as single or multiple fiber, and 
20 shown in Figure 1 5, along with a map of the energy in the odd channels which 
shows little structure. The magnitude of the single fiber channel, and the 
estimated maps of 0and <f>, is shown in Figure 14 for a different slice from the 
same data set. 

Discussion and Conclusion 

25 The recognition of the MR sensitivity to diffusion anisotropy [4] 

occurred shortly after the initial studies of isotropic diffusion [2, 3], but really 
blossomed with the recognition of its utility for the study of fibrous biological 
systems [10, 1 1], The natural first step in quantifying anisotropic diffusion is 
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to assume a Gaussian model for a single fiber, from which a diffusion tensor 
model ensues [8], From this model can be estimated the fiber diffusivities and 
orientation, and the conditions sufficient to capture this information [9]. While 
it was recognized that this method is restricted to single fiber [8, 12], 
5 development of a technique to measure more complex systems was only 
recently proposed [13]. With the general goal being to investigate complex 
anatomical structures with no a priori assumptions about the diffusion 
characteristics, the sampling criterion is no longer based upon a presumed 
model, and so is phrased in terms of being "unbiased" in the sense that no 

10 direction is assumed preferable. The result is sampling on a sphere with a 

radius defined by the 6-value with sufficient density to detect diffusion changes 
in different directions [13]. The question then arises: How does one 
characterize the diffusion from such measurements? One approach, suggested 
by Wedeen [14], is to extend the measurements to several radii by collecting 

15 "shells" of high angular spherical resolution data and Fourier transforming the 
data to produce a "q-space" image [26]. This method is sensitive to restricted 
diffusion because it samples a range of 6-values, but is unnecessarily 
complicated for the determination of non-Gaussian diffusion arising, for 
instance, from multiple fiber directions. As pointed out inn [15], the magnitude 

20 of b, in conjunction with high angular resolution sampling, is sufficient for this 
purpose. The more pertinent question is how to characterize HARD 
measurements made at high 6-vaIues, which is the subject of the present paper. 

Our approach to characterizing the diffusion measured with the HARD 
technique is based on characterizing the shape of the measured apparent 
25 diffusion coefficient. Using the methods of group theory, this characterization 
is that HARD measurements can be decomposed into irreducible 
representations of the rotation group in which isotropic, single fiber, and 
multiple fiber components are three, separable direct sum subspaces. In our 
initial investigation of this problem, deviations from a spherical surface in the 
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form of the variance of the measurements were used as a measure of 
anisotropy, the spherical diffusion variance [15]. While this had the advantage 
over the diffusion tensor method in identifying regions of anisotropy not well 
characterized by a single diffusion tensor, it is unable to quantify this 
5 anisotropy in any meaningful way. In particular, it does not allow the 

quantitation of either the magnitude or direction of diffusion. The present work 
can be seen as a formalization of this approach, since the spherical harmonics 
higher than 7 0 °, which contribute to the spherical diffusion variance, 
characterize how the various anisotropic components contribute to the variance. 

10 A strength of the approach is that it does not require any a priori 

information about the diffusion. The utility of the decomposition results from 
the group algebra imposed by the symmetries of both the measurement scheme 
and the diffusion. The decomposition allows distinction of diffusive and non- 
diffusive signal variations, as well as distinction amongst diffusion variations. 

15 In particular, the diffusion channels can be broken down into direct sum 

subspaces representing isotropic, single fiber, and multiple fiber components. 
Asymmetries produced by experimental artifacts fall into channels impossible 
to reach by diffusion, thereby providing a direct means of noise reduction 
within the diffusion channels as well as a means for identifying artifactual 

20 effects. 

The numerical computation of the spherical harmonic transform (SHT) 
was implemented in the most naive fashion, by direct computation of the 
coefficients by a discretized version of the integral Eqn 21, in a fashion 
analogous to a discrete Fourier transform (DFT) on a sphere. This method is 
25 time consuming. Unfortunately, there is no matrix decomposition of the SHT 
analogous to that used to implement the fast Fourier transform (FFT) from the 
DFT that would facilitate a fast SHT. Several methods for fast SHT 
computation have been proposed (e.g., [27, 28, 29, 30, 23, 24]) and will be 
investigated for application to the current problem. 
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A very general problem that arises is the determination of the voxel 
composition from the SHD. The solution to this problem hinges on the ability 
to estimate the parameters of a voxel, in particular, the number of fibers, their 
volume fractions, anisotropics, and orientations. Even in the simple example 
5 discussed above of two fibers, it was shown that it is not possible to determine 
all of these uniquely, as the volume fraction and relative fiber orientation both 
affect the higher order SHD components in a similar fashion. In practice, fiber 
configurations within a voxel may be much more complicated than the simple 
two fiber model, making the problem of parameter determination exceedingly 
10 complicated. 

The power of the SHD, however, is the fact that the identification of the 
existence of multiple fibers is not dependent upon making this distinction: 
multiple fibers of any sort show energy in the higher channels. Moreover, even 
in lieu of a particular fiber model, the SHD allows the shape of the diffusion 

1 5 measurements to be quantified by the coefficients of the SHT. These may be 
then used to reconstruct the diffusion structure in each voxel depending upon 
the user model. With constraints, such as on the number of fibers, the 
computational complexity can be reduced. An example is shown in Figure 16 
in which two fibers were assumed. One solution for the composition in both 

20 angles and volume fractions was estimated simply by trial and error in the 

simulation. However, one can imagine formalizing this process by choosing a 
limited maximum number of fibers and searching over number of fibers, the 
relative angular displacements between the fibers and the fiber volume 
fractions in order to estimate these quantities. Incorporation of other imaging 

25 information may augment this estimation. 

For the purpose of determining the significance of the power in the 
multiple fiber channels in order to determine whether a voxel actually contains 
multiple fiber directions, some comparison with the single fiber channel is 
required. For this purpose, a simple statistic for determining the significance of 
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the multifiber channel, the Fractional Multifiber Index, was introduced. While 
this is a natural measure, there is no indication of its optimality. A more formal 
probabilistic analysis needs to be undertaken to determine a method for 
determining the significance of energy amongst the channels. 

5 We remind here that what we have loosely referred to as " fibers" are 

really fiber bundles, within which the water movement that produces the 
diffusion signal is most likely complex. We have implicitly assumed that all 
the fibers that make up a bundle are essentially identical. Moreover, we have 
also made the assumption throughout this paper that there is no exchange 
10 between fibers, so that the signals from the individual fibers add independently. 
This will not be true in general Both of these assumptions are made here, as 
they are in most of the DTI literature, not for lack of recognition, but because 
the true nature of the diffusion signal can be exceedingly complicated and is 
beyond the scope of the current paper. 

15 The characterization of fibers in terms of even orders of L is a 

consequence of the symmetry. Experimental artifacts, such as eddy current 
effects, that do not possess such symmetry, are not confined to the even 
channels and therefore appear in the odd L channels. Energy in these channels 
therefore indicates the presence of non-diffusion effects. The SHD, by 

20 automatically separating out some fraction of the non-diffusion energy, 

effectively reduces the noise in the diffusion channels. This has great potential 
for use on systems for which artifacts are present. 

In practice, the issue of multiple fiber directions within a voxel is 
intimately tied to the image resolution: For higher resolution there will be 
25 fewer voxels with multiple directions. Nevertheless, there will always remain 
heterodirectional voxels at any resolution. Moreover, the penalty in SNR per 
unit time can make high-resolution diffusion imaging over a large region 
prohibitive. A post-processing scheme capable of accurately identifying and 
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quantifying multiple fiber voxels may lead to more efficient acquisition protocols. 

One important application of our method is its incorporation of multifiber voxels 
into fiber tract mapping schemes. This will require utilizing estimates of the individual 
5 fiber orientations and volume fractions determined from the SHD of individual voxels. 
Additional machinery to keep track of multiple possible pathways of fibers passing 
through such voxels will then be necessary. 

It is worth reiterating that the described method reduces to the standard diffusion 

1 0 tensor method in the presence of single fiber voxels, so no penalty of information is 

imposed by its usage. Rather, deviations from the DTI model due to artifacts or multiple 
fiber directions are readily extracted and quantified, allowing a more complete 
description of complex diffusion processes in tissues. More specifically, in one 
implementation of the above spherical harmonic decomposition (SHD), a method for 

1 5 analyzing magnetic resonance imaging (MRI) data may include representing MRI 
diffusion data of a body part as a summation of spherical harmonic functions and 
separating terms of the spherical harmonic functions to represent different diffusion 
effects including information on anisotropy in diffusion embedded in the MRI data. This 
separation of the spherical harmonic functions suppresses noise in extracted diffusion 

20 data contributed from non-diffusion effects. 

FIG. 17 further shows one implementation of the above SHD processing. MRI 
images are constructed from high angular resolution diffusion data. The gradient 
directions are used to determine measurement angles and spherical Voronoi areas in the 
MRI images. Next, the spherical harmonic functions at specified gradient angles are 

25 computed and integration measures are also computed from the spherical Voronoi areas. 
Subsequently, a spherical harmonic transform for each voxel is performed to obtain 
coefficients for the spherical harmonic functions; The coefficients of the spherical 
harmonic transform are used to extract information on diffusion anisotropy. 

Notably, the SHD provides a general method for characterizing diffusion in MRI. 

30 In implementation, it is desirable to have a simple quantitative measure of anisotropy that 
can be used to easily compare experimental results, such as to summarize the diffusion 
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anisotropy (DA) in regions of interest, and amongst subjects. In DTI, the most common 
such measures at the relative anisotropy index (RA) and the fractional anisotropy index 
(FA) which are simple scalar measures of anisotropy calculated from the estimated 
diffusion tensor. See, e.g., Basser, P.J., Pierpaoli, C, Microstructural and physiological 
5 features of tissues elucidated by quantiative diffusion tensor MRI, Journal of magnetic 
Resonance B, vol 111, pages 209-219(1996). These measures can be generalized to the 
more general tensor expansion of the SHD as follows. The RA is defined as the ratio 
RA = -yjD A • D A I -yjDj • Dj , where D A is the anisotropic component of D and Dj is the 

isotropic component. The FA is defined as the ratio FA = yjD A • D A / -JD D . Using the 

10 definitions of RA and FA, but with the more general SHD expansion of D, we can write a 
more general forms the general RA or GRA and the generalized FA or GFA: 



where the a lm are the coefficients in the SHD. These scalar measures easily make the 
additional information obtained form the SHD apparent in the increased visibility of 
white matter (WM) structures in areas of complex fiber orientation, such as where the 

20 corpus callosum crosses hemispheres. Both the GRA and the GFA provides a simple and 
useful metric to quantitate DA from the results of the SHD. Since the FA is the most 
prevalent index used in the DTI literature, we will use it as the metric for by which DA 
will be measured in the DTI, and will use its generalization to the GFA as the index for 
the SHD. The HARD data are a superset of the data required for DTI, so that the same 

25 data set can be analyzed using both DTI and HARD methods. 
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While the foregoing specification has been described with regard to certain 
preferred embodiments, and many details have been set forth for the purpose of 
illustration, it will be apparent to those skilled in the art that the invention may be subject 
to various modifications and additional embodiments, and that certain of the details 
described herein can be varied considerably without departing from the spirit and scope 
of the invention. Such modifications, equivalent variations and additional embodiments 
also intended to fall within the scope of the appended claims. 
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